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Abstract 

The objective of this challenge is to develop a data-based probabilistic model of uncertainty 
to predict the behavior of subsystems (payloads) by themselves and while coupled to a primary 
(target) system. Although this type of analysis is routinely performed and representative of 
issues faced in real-world system design and integration, there are still several key technical 
challenges that must be addressed when analyzing uncertain interconnected systems. For 
example, one key technical challenge is related to the fact that there is limited data on target 
configurations. Moreover, it is typical to have multiple data sets from experiments conducted at 
the subsystem level, but often samples sizes are not sufficient to compute high confidence 
statistics. In this challenge problem additional constraints are placed as ground rules for the 
participants. One such rule is that mathematical models of the subsystem are limited to linear 
approximations of the nonlinear physics of the problem at hand. Also, participants are 
constrained to use these models and the multiple data sets to make predictions about the target 
system response under completely different input conditions. 

Our approach involved initially the screening of several different methods. Three of the ones 
considered are presented herein. The first one is based on the transformation of the modal data to 
an orthogonal space where the mean and covariance of the data are matched by the model. The 
other two approaches worked solutions in physical space where the uncertain parameter set is 
made of masses, stiffnessess and damping coefficients; one matches confidence intervals of low 
order moments of the statistics via optimization while the second one uses a Kernel density 
estimation approach. The paper will touch on all the approaches, lessons learned, validation 
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metrics and their comparison, data quantity restriction, and assumptions/limitations of each 
approach. 

Keywords: Probabilistic modeling, model validation, uncertainty quantification, kernel density 

1. Introduction 

To address the Sandia Challenge problem, NASA Langley assembled a group of researchers 
to study the tasking document and to interpret the information provided. Initially, after 
reviewing the documentation, it was apparent that no single approach appeared to be more 
suitable to tackle the challenge problem nor did a solution exist in the open literature. Instead, 
several different solution approaches were sought. In spite of the tremendous efforts of the 
Sandia group in crafting the challenge document, questions regarding the appropriateness of 
metrics, and what specifically were the “no modeling” ground rules still persisted and prompted 
many internal discussions. The two main modeling questions were: what was an acceptable level 
of modeling, and was the use of physical parameters, e.g., mass, stiffness, and damping for 
uncertainty modeling permitted. In the end, the work that is reported here maintains the spirit of 
the challenge problem, in our view, with Langley’s interpretation of the Sandia task, which is 
provided explicitly later in the document. 

From the onset it was clear that the key to obtaining a solution to the problem was to build 
probabilistic models for the parameters of the subsystem that could later be sampled in order to 
study the target system. Techniques to accomplish this would lead to models that “match” the 
experimental data provided. Without diminishing the work of statisticians, it was felt that a new 
solution approach would have to be tailored for this problem. Details on the proposed methods 
are presented next. 

2. Problem Statement 
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In the interest of time, the reader is referred to the Sandia Challenge document in Ref. [1] for 
a detailed description of the problem formulation, goals, and objectives. However, because of the 
several possible interpretations of the document requirements, it is important to restate the 
requirements as interpreted by the Langley team. Rather than rephrasing our interpretation of the 
requirements, it is best to postulate a series of questions that this work will attempt to address: 1) 
Is a linear time invariant (LTI) representation of the 3 -degrees-of- freedom (DOF) subsystem 
adequate to describe the nonlinear response data, i.e., does this simple model adequately capture 
the physics of the problem? 2) What metric(s) properly quantifies the agreement between a 
subsystem model and the experimental data? 3) How are the metrics computed and what are the 
associated confidence levels? 4) What are the acceptance/rejection criteria of the model and is 
the model found to be adequate by the criteria? 5) How can one determine that a model based on 
calibration data is suitable for validation? 6) Is updating of the calibration model recommended 
before using it to evaluate the accreditation data, and if so what is the updating process? 7) What 
is the impact of having additional number of tests in the database? How is this benefit/impact 
quantified? 

Answering each of these questions in a comprehensive manner is impossible in the time 
allotted to work this problem. Nonetheless it is recognized that the main purpose of the 
challenge is to address each question. Consequently, as a compromise, each solution approach 
was evaluated on its ability to address each question. Many approaches were abandoned due to 
their inability to answer a particular aspect of the problem. 

As a final note, it is important for the reader to remember that the same definitions in the 
Sandia document are used in the paper. For example, the subsystem is a 3 mass/spring/damper 
system with inputs applied at mass 1 for the calibration data and mass 2 for the validation data. 
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This subsystem is referred to as the 3 degrees-of-freedom linear time invariant 3-DOF LTI. 

When this 3-DOF LTI system is attached onto a flexible beam, it is referred to as the 
accreditation system but when attached onto a higher fidelity beam model it is the referred to as 
the target system. 

3. Evaluation of the LTI Model Scripts Using Metrics 

In this section, the adequacy of the LTI representation of the 3-DOF subsystem is studied. 

Unlike problems where engineers start with the physics, formulate the problem, and then solve it, 
for the challenge problem a set of MATLAB scripts were provided to approximate the physics of 
the 3-DOF nonlinear system. The goal here is to assess the adequacy of the scripts with respect 
to the nonlinear response data. For this, a metric that depends on the approximated lower and 
upper bounds of the support of the maximum and minimum singular values was used. If H M is 
the Frequency Response Function (FRF) matrix of the model and //? is the one extracted from 
the experimental data, model adequacy requires 

Max <t( H m ( jco , d , )) > Max a(H T ( jco , )) (3.1.1) 

\/d i Vdj 

and 

Min <j(H m (jco, d l )) < Min a(H T (jco, djj) (3.1 .2) 

Vd[ Vdj 

where d, is the vector of uncertain parameters, co is the frequency, cr is the maximum singular 
value of the input-output matrix, and a is the minimum singular value. This metric will be used 
to evaluate the script for the 3-DOF LTI model provided. 

Although 60 calibration data sets were provided only 20 are statistically independent. For 
each experiment, the corresponding parameter data set was used with the script, i.e. H M , to 
simulate the 3-DOF system response and to compute Eqs.(3.1.1)-(3.1.2). Figure La. shows with 
a solid blue line the test data and with a dotted red line are results using the Sandia script for the 
calibration experiments, a for test and model are the top two lines and a_ are the bottom two 
lines. Discrepancies in the lower singular value bound prompted a close examination of the 
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Sandia scripts that revealed numerical problems with a matrix inversion operation. After 
replacing a numerical matrix inversion with its closed-form solution, the curves in Fig. l.b. were 
computed. 

Despite the fact that results with the modified script proved to be better, when comparing 
singular values, it is still difficult to judge how well the 3-DOF LTI model predicts time 
observations. For this another metric was defined in terms of the 95% confidence intervals of 
the mean value of the maximum acceleration. Because of the limited number of experimental 
data sets, a /-distribution was used to determine confidence intervals for both test and analysis. 
For adequacy assessments, whenever possible the number of data sets for test and analysis was 
kept the same. This metric is given by the bounds of the 95% confidence intervals of the 
maximum acceleration based on the /-distribution. This implies that 

P a -t s/yfn < a <a + 1 s/yfn =0.95 (3.1.3) 

where a is the mean of the maximum acceleration, / is the /-distribution scaling for a population 
size of n, and S is the sampled standard deviation. We use this metric to compare the models 
corresponding to the Sandia script and the modified script. Figure 2a shows the mean of the 
maximum acceleration from test on the abscissa and results for the original Sandia script along 
the ordinate for subsystem masses 1, 2, and 3, as labeled. Superimposed at each estimated mean 
value are the 95 % confidence intervals plotted as thin bars for both test and the script; horizontal 
bars for tests and vertical bars for analysis. Data points for the mean of the maximum 
acceleration on the diagonal line indicate exact matching between the test and the model. Same 
information for the modified script is shown in Figs. 2a and 2b. To contrast this with results 
shown earlier, the Sandia and the modified scripts yield similar results. Consequently, the need 
to use different metrics to assess adequacy of LTI models even in those cases where the model 
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predicts well a particular observation is extremely important. Hereafter, the modified script will 
only be used. 

4. Definition of Metrics Used for Evaluation 

In the Sandia challenge document, a certification metric was provided to assess the 

survivability of a subsystem under loads and uncertainty. However, to evaluate the matching 
between the experimental data and the subsystem model two additional metrics have been 
chosen; a frequency domain metric using singular values of the FRF similar to those used in 
Refs. [2]-[3], and one that compares maximum acceleration and the corresponding confidence 
intervals. To define Metric 1, let a(H) be a diagonal matrix of singular values of the FRF. The 
average singular value over n observations is 

Metric l = ~y]tr(<j(H(j(o,d l ))) (4.1.1) 

n /=1 

where tr (- ) is the trace operator. This metric is computed not only for the model’s 

FRF, H M (j co,d,) but also for the test H T (jco,d,). Since the singular values of the FRF’s are 

random processes dependent on the frequency, several of their moments can be calculated and 
used for comparison to ascertain model “correctness”. Another important aspect of this metric 
is the fact it is independent of the input form, e.g., broadband, short pulse, sinusoidal, etc. 
Consequently, it is focused on the LTI structure and the numerical solution approach and not on 
how a particular input propagates through the model. 

Although the singular value metric works well to study matching in frequency domain, it is 
also important to establish a time domain metric to get a more direct look at the certification 
criteria. For this, the metric introduced in Section 3, referred hereafter as Metric 2, will be used. 
Hence, Metric 2 is given by the bounds of the 95% confidence intervals of the maximum 
acceleration. Hence, 
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Metric 2 = 


(4.1.2) 


a —tS/ \fn,a +tS/ yfn 

5. Probabilistic Modeling 

In the following sections three probabilistic modeling methods are presented to create new 

parameter populations from experimental data observations. These are the methods used to 
address the certification problem later in the paper. 

5.1. Probabilistic Modeling Approach (PM-1) 

Although there are many aspects of the challenge problem that have to be carefully 
addressed, the most important aspect is the approach used to generate parameter values outside 
the given data set. Any solution approach must not only provide parameter data values beyond 
what was provided but more importantly these new parameter values should exhibit similar 
statistics. Although matching of the statistical features is a reasonable approach to follow, there 
are many dimensions of the matching process that could be examined. In this approach, the 
mean and covariance of the modal data are matched by a correlated normal random vector. This 
approach transforms the data from the parameter space to an orthogonal space using singular 
value decomposition. It is in this orthogonal space that a normal distribution is built. The 
correlation in the data is captured by designing the filter T. The steps to design the filter are 
explained next. 

First, the three data sets for high, medium, and low inputs are averaged in order to obtain 
an independent sample. To avoid numerical conditioning problems that often arise when 
parameters have drastically different magnitudes, the data are then normalized using the 
maximum value of the parameters. The normalized data will be referred to as the parameter 
matrix D, where D = [d l d 2 d n ] . The next step is to decompose the parameter matrix into 
principal components using singular value decomposition (SVD). Given a kxn parameter 
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matrix assembled from n observations of k parameters in the vector d, its SVD decomposition is 
written as; 

D = USr=U[s 0]M = t/ s 2 (5.1.1) 

kxn kxk kxn nxn kxk £ xw 

Since the number of observations n is greater than the number of parameters k, there are only k 
nonzero singular values in the matrix S and they are stored in the matrix partition s. 

Furthermore, Eq.(5.1.1) shows the partitioning of V to separate those columns associated with 
the nonzero singular values. Since the vectors in U form a k th dimensional orthogonal set, by 
using the transformation T = Us the problem can be transformed from the input parameter space 
d to an orthonormal space 77 using 77 . = T l c/ / , where // , is the / h column of 77 . In this 

transformed space it is possible to truncate the parameter space according to the singular value 
magnitudes. Although truncation is not performed in PM-1, both PM-2 and PM-3 employ 
reduced dimension uncertainty models. 

To generate a new population to approximate the random vectors in the columns of the 
transformed space 77 , define a new k x 1 random vector n e = Ty where y is a k x 1 random 

vector from a normal distribution with E\y \ = T x fj, , E[yy T ] = / , and T is a filter matrix yet to be 
defined. In this definition // is the expected value of the transformed parameter vector 
ju = E[r/j] and the data covariance in the transformed space is Q= E\ij i)J\ = rf 7 . This last 
equality is a result of a Cholesky decomposition of the data covariance matrix. Hence, using 
these definitions E\ n e ] = // , E[n e n T e ] = Q , and E[d\ = T/u . This process is similar to that of 

defining a fully correlated multi-variable normal distribution, but the added flexibility is in the 
ability to truncate the parameter space using singular values. 

5.2. Probabilistic Modeling Approach (PM-2) 
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This approach consists of two main steps, namely, the uncertain space truncation and the 
sizing of a joint Probability Density Function (PDF). In the former step, the dimension of the 
uncertain parameter space is reduced by approximating the dependencies among the uncertain 
parameters with linear combinations. This step is based on the SVD of the matrix of 
observations. Joint PDF sizing is carried out in the resulting space, by solving an optimization 
problem that minimizes the offset, to be defined later, between the observations and the 
predictions for the subsystem. Inequality constraints are used to prevent identifying models 
leading to infeasible subsystems. Details on the implementation of these steps are presented 
next. 

5.2.1, Truncation of the Uncertain Parameter Space 

In contrast to PM- 1 , which works with modal parameters directly, this approach uses 

physical parameters. Consequently, the row dimension of the parameter matrix is k=9 and not 
A=15 as presented in PM-1. The approach developed here will define the uncertainty space as a 
general r-dimensional space, with specific numerical values for r given later in the paper. As 
before, the high-, medium- and low-level excitation results for the same realization of the 
subsystem are first averaged to make a statistically independent sample. In contrast to the PM-1, 
this method normalizes the data using the estimated mean value of the parameters. The 
normalized data will be collected in the observation matrix X. As a result, the i th row vector of Xy 
contains the mean-normalized observations of the z th parameter. The SVD of this matrix leads to 

X = USV (5.1.2) 

SH Sr 1 S - 1 S - 1 x 7 

kxn kxk kxn nxn 

An r-dimensional approximation to this matrix is given by 

X = USV (5.1.3) 

V Sh S- 

lcxn kxr rxr rxn 
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where the underlined matrices result from extracting sub-matrices from Eq. (5.1.2) in the referred 
positions, e.g. U has the same first r columns of U. This approximation is equivalent to setting 
the smallest non-zero k-r singular values of the matrix S in Eq.(5. 1 .2) equal to zero. The 
approximation error due to truncation, given by 

is used in the truncation process to select r. 

Note that Eq. (5.1.3) contains a linear combination of r uncorrelated parameters (the ones in 
V), possibly dependent among themselves, that approximates the k parameters of interest (the 
ones in X). Hereafter, we will refer to the £ parameters in the left hand side of Eq.(5.1.2) as the 
Target Parameters and to the r parameters of V in Eq. (5.1 .3) as the Slack Parameters . While the 
target parameters might have a clear physical interpretation, e.g. they could be in either the 
modal or the physical space, the slack parameters do not. 

If an r-dimensional probabilistic model of the slack parameters is available Eq.(5.1.3) can be 
used to generate samples of the target parameters. The resulting subsystem model however, must 
be able to prevent the occurrence of infeasible subsystem realizations. This requirement is 
achieved by parameterizing the slack parameters in terms of r target parameters. Such a step will 
be referred to as the re-parametrization. Its implementation is as follows. From Eq. (5.1.3) select 
an arbitrary set of r row equations out of the k available. This set can be written as 

B = C S V (5.1.5) 

u '-' ^ s= ^ v ’ 

rxn rxr rxr rx „ 

where B is extracted from X and C is extracted from U and S. If C is a non-singular matrix 
Eq.(5.1.5) provides the desired re -parameterization. Combining Eqs. (5.1.2) and (5.1.5) one can 
obtain the approximation 
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(5.1.6) 


X = U C 1 B 

^ ‘ -v_l ^ 

kxn rxr rxr rxn 

This expression provides a parameterization of the k target parameters in terms of an r- 
dimensional subset. Since the r parameters in B have a clear physical interpretation, one can 
easily construct a sound probabilistic model for them. 

Once r rows of X are selected to form B, the independent parameters in Eq.(5. 1 .6) are 
prescribed. If such parameters are not a basis for the corresponding space, approximation error is 
introduced. In practice, this error will always be present. In order to minimize such an error, we 
select the combination of r target parameters that makes the condition number of C as close as 
possible to one. The closer the condition number is to one the smaller the approximation error. 

5.2.2, Joint PDF Sizing 

Equation (5.1.2) and a //-dimensional probabilistic model of uncertainty for the parameters in 
V, constitute one type of uncertainty model for the subsystem paramters. Similarly, Eq. (5.1.3) 
and a r-dimensional probabilistic model of uncertainty for V constitute another type, one for 
which truncation error is present. Another type results from using Eq. (5.1.6) and a r-dimensional 
probabilistic model of uncertainty for B. The resulting model for this type has both truncation 
and re -parameterization error. 

In this section we will consider the model based in Eq. (5.1.6). Hence, the only task that 
remains to be done is to determine a r-dimensional joint PDF for the uncorrelated, possibly 
dependent, parameters in B. Assume that a sample of the r parameters in B, namely b, is given 
by 

b=Ah(q) (5.1.7) 


where A is a transformation matrix, and h(q ) is a sample of an r-dimensional vector of 
independent random variables. The vector q is used to explicitly refer to the parameters that 
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specify the PDF of h while the matrix A is used to allow for correlation of the independent 
variables. However, the free parameters A, q and the structure of h are to be specified. 

If the target parameters are physical parameters, the components of b must be strictly 
positive. This requirement is attained by selecting a matrix A whose components are non- 
negative, and by choosing a PDF for h whose support can only assume non-negative values. 

Once the structure of h is prescribed, e.g. h is lognormal, we have r unknowns and r~ 
constraints from A, and dim(q ) unknowns from h. In what follows, we assume that h is a vector 
of independent lognormal random variables. Hence, dim(q)=2r. 

If Eq. (5.1.6) is rearranged such that the first r target parameters in X are the ones in B ; joint 
PDF sizing is performed by solving the following optimization problem 

< A,q>= axgmm{p T Wp : A ij > 0 for i = l,...,rand j = 1 

UC~ 1 A i/ >Ofori = r + l,...,kandj = l,...r} 

where p is a vector composed of metrics that quantify the offset (to be introduced later) between 
the observations and the predictions of the model, and W is an arbitrary, strictly positive, 
weighting matrix. The constraints in Eq. (5.1.8) are solely used to eliminate models leading to 
physically unrealizable subsystems. 

There are many ways to quantify the “offset” in p. The metrics considered are introduced 
next. In what follows we will denote by x 0 the /'-dimensional vector of target parameters in B 

corresponding to the observations and by x m the one corresponding to the model. Note only r 
target parameters are used to size the joint PDF. 

5.2.3. Low Order Statistics 

Arbitrary order statistical moments of the target parameters can be readily approximated from 
the observations. However, since the population size is quite small, only the first few moments 
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are meaningful. For this reason we have only considered the first, second and third order 
moments of the target parameters. On the other hand, the predicted moments resulting from 
using Eqs.(5.1.6)-(5.1.8) and any given q, can be calculated analytically. Details on their 
calculation are omitted due to space limitations. The metric of interest is 


offset ■ = wf 


IIQ[x 0 ]-E[x m ]|| 

||Q[x 0 x 0 t ]-E[x m x m t ]|| 


(5.1.9) 


ZllQ[ X Oj X 0 Tx «]- E [ X Mj X M Tx M]ll 

_ J =1 

where wi > 0 is a vector of weights, Q[-] is the sampled-mean operator, E[-] is the expected value 
operator, and ||-|| is the Euclidean norm. Note that this offset does not take into account the 
population size, e.g. the fact that for the calibration experiment there are only n = 20 statically 
independent observations available. 

5.2.4. Confidence Intervals of the Low Order Statistics 

This metric quantifies the offset between the 95% confidence intervals of the statistics in 

Eq.(5.1.9). To implement this, n (the same number of observations available) Hammersley 
Sequence samples (HSS) [4] of the subsystem model are generated. This sampling method is 
used since the resulting samples not only characterize well the joint PDF but also because they 
are deterministic. Note that these confidence intervals measure the sampling PDF of the 
statistics, which is assumed to be Gaussian. The Maximum Likelihood Estimator method (MLE) 
is used to calculate these intervals. 

If the 95% confidence interval of the random variable a is denoted as [f(a) ,y + (a)], the 
corresponding metric for this offset is 
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offset , = wf 


llr(Q[xo])-r(Q[x M ])IHIr + (Q[xo])-r + (Q[x M ])ll 

ll 7 '(Q[x 0 x 0 T ])-f‘(Q[x M x M T ])||+||^ + (Q[x 0 x 0 T ])-f + (Q[x M x M T ])|| 

Zll 7 '(Q[x 0j x 0 x 0 T ])-/-(Q[x Mj x M x M T ])||+||^ + (Q[x 0j x 0 x 0 T ])-^ + (Q[x Mj x M x M T ])|| 


(5.1.10) 


where w? >0 is a vector of weights. If more observations are available, smaller confidence 
intervals for the statistics result. This will reduce the freedom available for sizing the joint PDF 
that exist when n is small, i.e. the offsets corresponding to different subsystem’s models differ 
more as the number of observations increases. 

5.3. Probabilistic Modeling Approach 3 (PM-3) 

Kernel Density Estimation (KDE) provides another technique to approximate unknown 
parameter distributions functions from observations. Due to space limitations, the technical 
details of the KDE approach, found in Ref. [6], will not be discussed here. Instead, potential 
users are referred to the MATLAB-based program called the Kernel Density Estimation Toolbox 
described in Ref. [7], which is a Multivariate KDE software package that was utilized in this 
paper. The general appeal of the KDE approach is that it does not impose a parametric form on 
the density function estimate prior to any actual estimation. 

Similar to PM-2, this approach operates on data in the physical domain (mass, stiffness and 
damping) as opposed to the modal domain (frequency, modal damping, and mode shapes) as in 
Probabilistic Model 1 (PM-1). Furthermore, analogous to PM-2, the KDE approach in PM-3 will 
employ the singular value decomposition presented in Eqs.(5.1.2) and (5.1.3) to reduce the 
dimensionality of the uncertainty space. For this, the error metric in Eq.(5.1.4) is also used to 
control truncation error when selecting r in Eq. (5. 1 .3). Again, as in PM-2, avoidance of non- 
physical parameter values must be guaranteed. In PM-3, this is achieved using a brute-force 
method of sampling the identified kernel density estimates using very large sample sets and 
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rejecting those density estimates that produced non-physical values. There has been a great deal 
of research on kernel density estimation and it is well-known that the selection of the actual basis 
functions is not as critical as is the selection of the basis function’s bandwidth, see Ref. [6]. 
Nonetheless, the KDE Toolbox does provide user-selectable basis functions as well as a variety 
of methods to prescribe their bandwidths. The results obtained in this work use Gaussian basis 
functions with bandwidths computed using the plug-in approach, see Ref. [7] toolbox, for 
minimizing an approximate mean integrated square error criterion. 

The steps used in generating the kernel density estimates area as follows: 

1 . The raw modal data is processed to provide physical parameters (M, C, and K). This is 
done separately for both calibration and validation data sets for all force levels. 

2. Physical data for the three force levels is averaged to provide one parameter set for each 
of the twenty system realizations, resulting in twenty data sets using calibration data and 
twenty data sets using validation data. 

3. Singular value decomposition can be performed on the aggregate [M,C,K] combined 
dataset generated in step 2, or separately for [M], [C], and [K], Regardless, the dimension 
of the uncertainty space is selected using the error in Eq. (5.1.4). Note that the data for 
step 2 may be treated separately, i.e., two sets of twenty points, one for calibration and 
one for validation, or as a single forty point dataset (calibration and validation combined). 
Both approaches are used. 

Once the data had been properly processed, the KDE is performed as outline above. 

6. Probabilistic Modeling Results 

In the following, results from the PM-1, PM-2, and PM-3 approaches are used to predict the 
behavior of the accreditation and target system, but first results from each approach with the 
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subsystem only are discussed. Also discussed is the evaluation of the certification criterion by 
sampling the resulting subsystem models. 

6.1. Probabilistic Method-1 (PM-1) 

Recall that this approach matches mean and variances of the experimental data. To give an 
idea of the type of matching achievable with this approach, Fig. 3a shows six plots, top three 
leftmost show mode 1, 2, and 3 frequencies, and bottom three are the corresponding damping 
values. Star symbols correspond to calibration and validation data while superimposed are lines 
for the 99.9% confidence intervals obtained using a maximum likelihood estimate. Figure 3b 
shows the experimental covariance matrix E[r/ / 77 ^ ] in transformed coordinate space at the top 

and the computed covariance using 40 data sets at the bottom. As expected, the computed 
covariance matrix does not match exactly that of the experimental data set for small sample sets. 
In order for the two to agree, an infinite number of samples would be required. 

A better way to evaluate a PM approach is to sample the resulting subsystem model and to 
evaluate the sample mean and representative percentiles of Metric 1 .Fig. 4 shows these results 
for a parameter population of 200 samples. In here, the solid blue line is the experimental 
average from 20 calibration data sets whereas in dashed red is the sampled mean of the model 
and dashed greens are the corresponding 5 and 95 percentiles, i.e. if F denotes a sampling-based 
approximation to the Cumulative Distribution Function (CDF) of the random process associated 
with the singular values, the percentiles l x (co) and l 2 (co) are curves satisfying F(l t ) = 0.05 and 
F(l 2 ) = 0.95 . Similarly Fig. 4b, shows results for the combined 20 calibration and 20 validation 

data sets; solid blue is the averaged for the experiments, dashed red is the model mean, and 
dashed green are the 5 and 95 percentiles. A subtle difference in the computation from Fig. 4a 
and 4b worth mentioning is the fact that in Fig. 4b, the singular value metric contains data from 
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two different inputs. A quick view of the plots shows that PM-1 and the 3-DOF LTI model 
predicts the average experimental results very well and that the percentiles are able to explain 
observations even in the presence of significant variations in the experimental response. 

6.2. Probabilistic Method-2 (PM-2) 

The approach presented in Section 5.2 was applied to the structural dynamics challenge, 
and the main observations and results are presented next. Since the structure of the subsystem is 
known, working in the physical parameter space allows one to reduce the dimension of the 
uncertain parameter space from 12 to 9, i.e. from C , and fa for z=l ,2,3 with the three mass 

normalized equality constraints for the eigenvectors, to m,, k t and c, for i= 1,2,3. The truncation 
error of Eq. (5.1.4)for the calibration and validation data and the random vibration experiments 
led us to chose r=4. At such value, e=0.0426 and the best set for re-parameterization is given by 
m 2 , C 2 , A'i and A 3 . 

Several joint PDFs were sized using different Ws in Eq.(5.1.8). Due to space limitations, we 
only present here results corresponding to one for which the offset metrics in Eqs.(5.1.9) and 
(5.1.10) were combined. 

In the figures that follow, information for all target parameters is displayed even tough the 
sizing of the joint PDF was performed using only four of them. In other words, while Eq. (5.1.8) 
only uses information on m 2 , C 2 , k\ and A 3 , the figures show statistics and confidence intervals 
corresponding to the 9 physical parameters. Calibration and the validation data were used to 
generate PM-2. Figure (5) shows the low order statistics for the observations (line with dots) and 
the model (line with circles). A small offset between the observations and the model is attained. 
Figure ( 6 ) shows the confidence intervals for the standard deviation of the statistics in Figure (5) 
using the same line convention. A good matching of the confidence intervals was also attained. 
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Note that a good matching in the sense of Eq.(5.1.9) does not imply the one of (5.1.10), and vice 


versa. 

As with PM-1, PM-2 is now compared to the experimental data using Metric 1. Figure 7a 
shows Metric 1 for calibration data and calibration/validation in shown in 7b for a parameter 
population of 200. As before, the solid blue line corresponds to the mean of the experimental 
data, dashed red is the model mean, and dashed green lines are the 5 and 95 percentiles. Similar 
to results shown in Fig. 4, the average from experiments is predicted well by PM-2 and the 
percentiles indicate that this model would also capture significant off-nominal experimental 
observations. 

6.3. Probabilistic Method 3 (PM-3) 

The KDE method, as presented in Section 5.3, was used to model the uncertain 3-DOF FTI 
subsystem data. The approach taken was to develop a KDE model for each of the physical 
parameters independently, i.e., treat [M], [C], and [K] separately. Employing Eq. (5.1.4) to 
determine the effect of truncation error results in a one-dimensional subspace for K and two- 
dimensional subspaces for M and C, i.e, [rk,r m ,r c ]=[ 1,2,2], for a total of five parameters. 
Comparisons of the normalized identified KDE models parameters versus those from the 
calibration data are shown in Fig. 8, with differences less than one percent. 

As before, the PM-3 model is evaluated using Metric 1. Figure 9 shows the corresponding 
results using the same convention of Figs. 4 and 7 for the calibration data. It is seen that the 
mean of KDE model tracks test mean (in blue) well, but more importantly stays within the 
percentiles throughout the frequency range of interest 

7. Discussion of results for target and accreditation system 

Until now the discussion focused on the 3-DOF LTI subsystem model and the development 

of probabilistic models for the parameters. However, it is important to remember that the goal of 
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the challenge problem is to make predictions, using what was learned from the subsystem to 
predict failures when integrated into the accreditation and the target system. Assessments for 
this, based on PM-1, PM-2, and PM-3 are presented next 

To begin studying the accreditation system, adequacy of the integrated model provided by 
Sandia must be addressed first, where the model now includes the 3-DOF LTI, PM-1, or PM-2, 
or PM-3, and a model of a flexible beam. Due to space limitations, the results in the following 
are from probabilistic models obtained using the calibration data only, but results with both 
calibration and validation are comparable. Figure 10a and b shows results for Metric 1 and 2, 
respectively, using PM-1 with a parameter population arbitrarily chosen to be 200. Results for 
Metric 1 are shown in Fig. 10a with the solid blue lines corresponding to data from three 
experiments, whereas dashed lines correspond to the mean (red) and the 5 and 95 percentiles for 
the model. It is reasonable to assume that if the model is to fully explain the test observations all 
test results must lie within the arbitrarily selected percentiles bounds. This is the case for most of 
the frequency points To contrast, time domain results for Metric 2 are shown in Fig. 10b with 
the abscissa showing the mean of the maximum acceleration for test versus model results plotted 
along the ordinate. To show uncertainty, lines from the mean, labeled ‘o’, shows the 95 % t- 
distribution confidence interval of the mean; horizontal lines for test and vertical lines for 
analysis. Since the analysis results are obtained with a population of 200, the uncertainty bounds 
are significantly smaller than test, which only has 3 data sets. Also note that the results for 
maximum acceleration of Mass 3 are close to the diagonal, indicating excellent matching of 
mean value predictions with tests. 

Results for Metrics 1 and 2 using PM-2 and PM-3 approaches are shown in Figs. 1 1 and 12. 
Amazingly, in spite of the drastically different approaches, when compared through the output of 
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the accreditation system the same conclusions can be drawn as for PM-1. However, both PM-2 
and PM-3 are able to realize a parameter population using a 4 th and a 5 th dimensional 
probabilistic model as opposed to 15 in PM-1, hence, significantly reducing the dimensionality 
of the uncertainty space. 

Finally, to investigate survivability of the subsystem as part of the target system, PM-1, PM- 
2, and PM-3 are all used along with a prescribed excitation signal to compute the failure 
probability as stated in the Sandia document: 

P f = P[max \a 3 (t)\ > 1.8x 10 4 m/sec 2 ] (7.1.1) 

The CDF of the maximum absolute value of the acceleration of mass three will be used to 
evaluate the models from calibration, validation, and accreditation data and to make a prediction 
on the target system response. Results for the case using Calibration data only are shown in 
Figure 13 in terms of the CDF for PM-1 in 13a, PM-2 in 13b, and PM-3 in 13c. With all three 
methods a parameter population of 5000 is used and outputs are generated for each case. As 
noted in the figure the failure probabilities predicted with PM-1 is 0. 17, PM-2 is 0.22, and PM-3 
is 0.20. But no significant difference in the shape of the CDF nor the in support set of the 
random variable is noticeable. Consequently, when the probability values are compared to the 
accreditation requirement that P f< 10 2 , all three methods predict failure of the subsystem on 

the target system. Figure 14 shows the CDF results for the case when both calibration and 
validation data are used. In contrast to the results in Fig. 13, the failure probabilities predicted 
with PM-1 is 0.13, PM-2 is 0.12, and PM-3 is 0.20. Although slightly different from the case 
with calibration data only, this case also predicts failure of the subsystem. Also noticeable is 
how similar the CDF are among themselves and as compared to the ones in Fig. 13 and 14. 
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However, the failure probability estimate corresponding to PM- 1 and PM-2 are considerably 
more sensitive to the inclusion of the validation data than the one from PM-3. 

8. Concluding Remarks 

Three approaches to parameter uncertainty modeling, based on experimental data, have been 
proposed and applied to the Sandia Challenge problem. To assess the adequacy of the model 
structure in the Sandia scripts and our ability to predict experimental observations, two metrics 
were defined and used to validate models at the subsystem and system levels. When using a 
frequency response metric a numerical deficiency of the Sandia script surfaced, but later it 
proved not to significantly impact the maximum acceleration. Note however, that this single 
metric does not imply a good matching in the time-domain per se. Amazingly, all three 
approaches produced almost identical CDFs for the acceleration performance metric and 
consistently showed that the regulatory requirement for the target system was violated. 
Nonetheless, using PM-2 and PM-3 the number of parameters needed to describe the uncertainty 
was reduced from 15 to 4 and 5. 

For this particular problem, all three probabilistic methods, created using calibration data 
only, provided consistent results and predicted failure of the target system. The inclusion of the 
validation data did not substantially change the assessments resulting from using calibration data 
only. Finally, the fact that all three methods produced similar results is only an indication that 
data quantity limitations did not allow for the proper exploitation of strengths and weaknesses of 
each technique. Clearly, PM-2 and PM-3 reduction in the number of parameters required to fully 
characterize the uncertainty is remarkable. 
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la) Sandia original scripts lb) Modified scripts 

Fig. 1 Comparison of calibration data singular values bounds using the original 
Sandia script and revised script. 
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2a) Sandia original scripts 2b) Modified scripts 

Fig. 2 Comparison of Metric 2 for the calibration data using original Sandia 
script and revised script. 
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3a) Sample frequencies and damping 


3b) Covariance 


Fig. 3 PM-1 results with 99.9 % MLE Cl Using Calibration and 
Validation Data 
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Fig. 4 PM-1 Comparison of Metric 1 (Population size 200) 
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Fig. 5 1st, 2nd, and 3rd moment results from matching PM-2 parameters with 
Calibration data 
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Fig. 6 Confidence intervals for the standard deviation of the results in figure 5 
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Metric 1 




7a) Calibration Data 


7b) Calibration & Validation 


Fig. 7 PM-2 Comparison of Metric 1 (Population size 200) 
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Fig. 8 Identified KDE parameters versus experimental data 
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9a) Calibration Data 9b) Calibration & Validation 

Fig. 9 PM-3 Comparison of Metric 1 (Population size 200) 




4.5 
4 

3.5 
3 

2.5 
2 

1.5 

Test Max Accel x 10 4 (in/s?) 


v 

•f - r 7* 

i i / 

— 

- Mass 3 . i 

— 

l ! J 

Mass 2 

i t 


i J. 

/ ' 

/ 1 1 

/ \ 1 

Mass 1 | 

: * 

1 1 

1 1 

1 1 


2 3 4 


10a) Metric 1 


10b) Mean and 95% Cl for Metric 2 


Fig. 10 Verification of Accreditation System Response Bounds PM-1 
(Population Size 200) 
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Metric 1 



11a) Metric 1 
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Fig. 1 1 Verification of accreditation system response bounds PM-2 
(Population Size 200) 



in 

£ 


4.5F 


3.5 


— 




Mass 3 ~ 

7 

+ 


1 

A 



Mass 2 

+ 



' 

-1- 


Ma 

ss 1 







) f 



2 3 4 

Test Max Accel x 10 4 (in/S 2 ) 


12a) Metric 1 


12b) Mean and 95% Cl for Metric 2 


Fig. 12 Verification of accreditation system response bounds PM-3 
(Population Size 200) 
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13a) PM-1, Pf=0.17 


13b) PM-2, P t =0.22 


13c) PM-3, P t =0.20 


Fig. 13 Target System Mass 3 maximum acceleration CDF using Calibration Data 
(Population Size 5000) 
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Fig. 14 Target System Mass 3 maximum acceleration CDF using Calibration and 
Validation Data (Population Size 5000) 
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